'''
Created on Jul 12, 2011

@author: oabalbin
'''

import sys
import numpy as np
from optparse import OptionParser
from collections import deque, defaultdict

def read_file(ifile, ofile):
    '''
    '''
    ifile = open(ifile)
    ofile =open(ofile,'w')
    
    for l in ifile:
        if l.startswith('#'):
            ofile.write(l)
            continue
        fields=l.strip('\n').replace(',','$').split('\t')
        DP4=fields[7].split(';')[3]
        DP4=DP4.replace('DP4=','').split('$')
        ct=0
        for i in DP4:
            ct+=int(i)
        if ct != 0:
        #if fields[4] != "." and fields[3]!=fields[4]:
            ofile.write(l)
    
    ifile.close()
    ofile.close()



def read_file_2(ifile, ofile):
    '''
    chr1    14824   .       A       .       33      .       DP=1;AF1=8.529e-05;CI95=1.5,0;DP4=1,0,0,0;MQ=23 PL
    '''
    ifile = open(ifile)
    ofile =open(ofile,'w')
    h=["SNV_ID","QUAL","VARIANT","DPREF","DPALT","MQ"]
    ofile.write(",".join(h).replace(',','\t')+'\n')
    
    for l in ifile:
        if l.startswith('#'):
            continue
        fields=l.strip('\n').replace(',','$').split('\t')
        chr,pos,qual,ref,alt=fields[0],fields[1],fields[5],fields[3],fields[4]
        info = fields[7].split(';')
        dpc = info[3].split('=')[1]
        ref_support, alt_support = int(dpc.split('$')[0]) + int(dpc.split('$')[1]),int(dpc.split('$')[2]) + int(dpc.split('$')[3])
        map_qual=info[4].split("=")[1]
        
        if alt==".":
            variant=str(0)
        else:
            variant=str(1)
        tl=[chr+'@'+pos,qual,variant,str(ref_support),str(alt_support),map_qual]
        ofile.write(",".join(tl).replace(',','\t')+'\n')
    
    ifile.close()
    ofile.close()

def read_file_3(ifile, ofile):
    '''
    chr1    881614  .       T       G       65.79   .       AC=0;AF=0.00;AN=2;DP=81;Dels=0.00;HRun=1;HaplotypeScore=6.81;MQ=58.60;MQ0=0     GT:AD:DP:GQ:PL  0/0:54,23:52:35.80:0,36,1506
    0/0:54,23:52:35.80:0,36,1506
    '''
    ifile = open(ifile)
    ofile =open(ofile,'w')
    h=["SNV_ID","QUAL","VARIANT","DPREF","DPALT"]
    ofile.write(",".join(h).replace(',','\t')+'\n')
    
    for l in ifile:
        if l.startswith('#'):
            continue
        fields=l.strip('\n').replace(',','$').split('\t')
        chr,pos,qual,ref,alt=fields[0],fields[1],fields[5],fields[3],fields[4]
        dest=fields[9].split(':')
        genotype1,genotype2=dest[0].split('/')[0],dest[0].split('/')[1]
        ref_support,alt_support=int(dest[1].split('$')[0]),int(dest[1].split('$')[1])
        
        if (alt_support)==0:
            continue
        if genotype2 =="0":
            variant=str(0)
        else:
            variant=str(1)
            
        tl=[chr+'@'+pos,qual,variant,str(ref_support), str(alt_support)]
        ofile.write(",".join(tl).replace(',','\t')+'\n')


def read_file_4(ifile, ofile):
    ''''
    #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SAMPLE
chr1    861305  chr1@861305     A       C       3.01    .       AA=M;SBALT=1;DPREF=37;DPALT=1;QBREF=35;DP=38;SBREF=2;MQALT=60;SBREF_Minus=11;SBALT_Plus=26;AF
=0.0263;MQ=58;SBALT_Minus=11;SB=2.17;SBREF_Plus=26;MQREF=58;QBALT=28    GT:AD:DP:GL:GQ  0/0:0,0:0:0,0,0:0
    
    '''
    ifile = open(ifile)
    ofile =open(ofile,'w')
    h=["SNV_ID","QUAL","VARIANT","DPREF","DPALT"]
    ofile.write(",".join(h).replace(',','\t')+'\n')
    
    for l in ifile:
        if l.startswith('#'):
            continue
        fields=l.strip('\n').replace(',','$').split('\t')
        chr,pos,qual,ref,alt=fields[0],fields[1],fields[5],fields[3],fields[4]
        info = fields[7].split(';')
        ref_support, alt_support = int(info[2].split('=')[1]), int(info[3].split('=')[1])
        
        if float(qual) <= 30.2:
            variant=str(0)
        else:
            variant=str(1)
        tl=[chr+'@'+pos,qual,variant, str(ref_support),str(alt_support)]
        ofile.write(",".join(tl).replace(',','\t')+'\n')



#ifile='/exds/users/oabalbin/projects/exomes/NA12878/analysis/sam_calls_quick/NA12878benign.sam_snps.mpileup.all.vcf'
ifile='/exds/users/oabalbin/projects/exomes/NA12878/analysis/sam_calls_quick/NA12878benign.sam_snps.mpileup.all.onlyvariants_isec_SS50Mh.vcf'
ofile='/exds/users/oabalbin/projects/exomes/NA12878/analysis/sam_calls_quick/NA12878benign.sam_snps.mpileup.all.onlyvariants_isec_SS50Mh.txt'
#read_file(ifile, ofile)
read_file_2(ifile, ofile)
#ifile='/exds/users/oabalbin/projects/exomes/NA12878/analysis/gatk_calls/NA12878.gatk_snps.b.raw.all_isec_SS30Mh.vcf'
#ofile='/exds/users/oabalbin/projects/exomes/NA12878/analysis/gatk_calls/NA12878.gatk_snps.b.raw.all_isec_SS30Mh_2.txt'
#read_file_3(ifile, ofile)
#ifile='/exds/users/oabalbin/projects/exomes/NA12878/analysis/varscan_calls/NA12878_snps_varscanV2.5_single_genome_07_01_2011_MQ20_isec_SS30Mh.vcf'
#ofile='/exds/users/oabalbin/projects/exomes/NA12878/analysis/varscan_calls/NA12878_snps_varscanV2.5_single_genome_07_01_2011_MQ20_isec_SS30Mh.txt'
#read_file_4(ifile, ofile)
